home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / test1.i < prev    next >
Text File  |  1995-07-26  |  14KB  |  399 lines
  1. /*
  2.    TEST1.I
  3.    Poorly vectorizing physics problem for timing Yorick.
  4.  
  5.    $Id: test1.i,v 1.1 1993/08/27 18:50:06 munro Exp $
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /* WARNING -- this code assumes 1-origin indexing is the default */
  11.  
  12. func test1(npass)
  13. /* DOCUMENT test1
  14.          or test1, npass
  15.      Track a mock "ablation front" as it propagates through a mesh.
  16.      If NPASS is given, the calculation is repeated that many
  17.      times.  The zoning, densities, temperatures, pressures, and
  18.      velocities are all computed arbitrarily, but the number of zones
  19.      and groups are taken to be representative of a typical 1-D
  20.      ablation calculation.  */
  21. {
  22.   if (is_void(npass) || npass<=0) npass= 1;
  23.  
  24.   now= split= array(0.0, 3);
  25.   timer, now;
  26.   n= npass;
  27.   while (n--) {
  28.     prevScale= prevAbl= 37;
  29.     for (i=1 ; i<=ntimes ; i++) accum, i;
  30.   }
  31.   timer, now, split;
  32.   timer_print, "Time per pass", split/npass, "Total time", split;
  33.  
  34.   if (prevAbl!=15 || prevScale!=8 ||
  35.       nallof(approx_eq(rorab,
  36.                [2.673480885e-03,2.566922597e-03,2.422313787e-03,
  37.             2.211817497e-03,1.870390096e-03])) ||
  38.       nallof(approx_eq(times,1.5)) ||
  39.       nallof(approx_eq(rordx,3.109107238e-03)) ||
  40.       nallof(approx_eq(rorsx,2.946068988e-03)) ||
  41.       nallof(approx_eq(rorhx,2.916081687e-03)) ||
  42.       nallof(approx_eq(dmax,3.034383698e-01)) ||
  43.       nallof(approx_eq(zs,5.973585920e-03)) ||
  44.       nallof(approx_eq(vel,-1.662706089e-01)) ||
  45.       nallof(approx_eq(acc,-3.148890163e-03)) ||
  46.       nallof(approx_eq(scl,2.061472601e-03)))
  47.     write, "***WARNING*** values returned by accum are not correct";
  48. }
  49.  
  50. func approx_eq(x, y)
  51. {
  52.   return (abs(x-y)/(abs(x+y)+1.e-35))<1.e-6;
  53. }
  54.  
  55. #if 0
  56. CPU seconds per pass:
  57.             IDL    Yorick     Basis       DAP
  58. HP730        -      0.41       8.21      8.15
  59. Solbourne  0.57     1.00       20.1      22.4  (1st trial)
  60. Solbourne  0.57     0.90       19.5      20.9  (2nd trial)
  61. Note:  Times on Solbourne fluctuate by ~10%
  62.  
  63. Run on forum (Solbourne) 8/Dec/92
  64.      forum[4] time /usr/local/lib/idl/bin.sunos.4.1.sun4/idl
  65.      IDL. Version 2.3.0 (sunos sparc).
  66.      Copyright 1989-1992, Research Systems, Inc.
  67.      All rights reserved.  Unauthorized reproduction prohibited.
  68.      Site: 1491.
  69.      Licensed for use by: LLNL - X Division
  70.      IDL> .rnew /home/miggle/munro/Yorick/include/test1
  71.      % Compiled module: PEAK.
  72.      % Compiled module: DERIV.
  73.      % Compiled module: ACCUM.
  74.      % Compiled module: APPROX_EQ.
  75.      % Compiled module: TEST1.
  76.      % Compiled module: INTERP.
  77.      % Compiled module: UNCP38.
  78.      % Compiled module: $MAIN$.
  79.      IDL> test1,1
  80.      IDL> exit
  81.      1.1u 0.3s 0:28 5% 0+1264k 4+0io 5pf+0w
  82.      forum[5] !!
  83.      time /usr/local/lib/idl/bin.sunos.4.1.sun4/idl
  84.      IDL. Version 2.3.0 (sunos sparc).
  85.      Copyright 1989-1992, Research Systems, Inc.
  86.      All rights reserved.  Unauthorized reproduction prohibited.
  87.      Site: 1491.
  88.      Licensed for use by: LLNL - X Division
  89.      IDL> .rnew /home/miggle/munro/Yorick/include/test1
  90.      % Compiled module: PEAK.
  91.      % Compiled module: DERIV.
  92.      % Compiled module: ACCUM.
  93.      % Compiled module: APPROX_EQ.
  94.      % Compiled module: TEST1.
  95.      % Compiled module: INTERP.
  96.      % Compiled module: UNCP38.
  97.      % Compiled module: $MAIN$.
  98.      IDL> test1,21
  99.      IDL> exit
  100.      12.5u 0.2s 0:30 42% 0+1432k 0+0io 0pf+0w
  101.  
  102.  
  103.      forum[6] yorick
  104.       Yorick ready.  For help type 'help'
  105.      > #include "/home/miggle/munro/Yorick/include/test1.i"
  106.      > test1
  107.             Timing Category     CPU sec  System sec    Wall sec
  108.               Time per pass       1.010       0.000       1.020
  109.              Total time       1.010       0.000       1.020
  110.       -----Total Elapsed Times-----       1.670       0.140      20.270
  111.      > test1,21
  112.             Timing Category     CPU sec  System sec    Wall sec
  113.               Time per pass       1.001       0.001       1.004
  114.              Total time      21.020       0.030      21.080
  115.       -----Total Elapsed Times-----      22.710       0.180      51.760
  116.      > test1
  117.             Timing Category     CPU sec  System sec    Wall sec
  118.               Time per pass       0.990       0.000       0.990
  119.              Total time       0.990       0.000       0.990
  120.       -----Total Elapsed Times-----      23.720       0.190     163.440
  121.      > test1
  122.             Timing Category     CPU sec  System sec    Wall sec
  123.               Time per pass       0.990       0.010       1.000
  124.              Total time       0.990       0.010       1.000
  125.       -----Total Elapsed Times-----      24.740       0.200     171.340
  126.      > quit
  127.      24.7u 0.2s 3:35 11% 0+1088k 18+0io 64pf+0w
  128.  
  129.  
  130.      forum[7] basis
  131.      Basis    (basis, Version 921125)
  132.      Run at 10:06:44 on 12/08/92 on the forum    machine, suffix 10470x
  133.      Initializing Basis System
  134.      Basis 7.0
  135.      Initializing EZCURVE/NCAR Graphics
  136.      Ezcurve/NCAR 2.0
  137.      Basis> echo=0
  138.      Basis> read /home/miggle/munro/Yorick/include/test1.bas
  139.      End of input file /home/miggle/munro/Yorick/include/test1.bas
  140.      Resuming input from TERMINAL
  141.      Basis> test1(1)
  142.  
  143.     CPU (sec)   SYS (sec)
  144.        20.150       0.000
  145.      Basis> end
  146.  
  147.     CPU (sec)   SYS (sec)
  148.        23.117       0.467
  149.      23.1u 0.4s 1:12 32% 0+2408k 20+0io 229pf+0w
  150.  
  151.  
  152.      forum[8] dap
  153.      DAP> read /home/miggle/munro/Yorick/include/test1.bas
  154.      DAP> test1(1)
  155.  
  156.     CPU (sec)   SYS (sec)
  157.        22.383       0.017
  158.      DAP> end
  159.      26.7u 1.1s 1:07 41% 0+2784k 17+0io 209pf+0w
  160.  
  161.  
  162. Run on tonto (HP730) 7/Dec/92
  163.      tonto[8] basis
  164.      Basis    (basis, Version 921125)
  165.      Run at 16:06:31 on 12/07/92 on the tonto    machine, suffix 2545x
  166.      Initializing Basis System
  167.      Basis 7.0
  168.      Initializing EZCURVE/NCAR Graphics
  169.      Ezcurve/NCAR 2.0
  170.      Basis> echo=0
  171.      Basis> read /home/miggle/munro/Yorick/include/test1.bas
  172.      End of input file /home/miggle/munro/Yorick/include/test1.bas
  173.      Resuming input from TERMINAL
  174.      Basis> test1(1)
  175.  
  176.     CPU (sec)   SYS (sec)
  177.         8.260        .010
  178.      Basis> test1(10)
  179.  
  180.     CPU (sec)   SYS (sec)
  181.        82.090        .020
  182.      91.7u 0.1s 4:07 37%
  183.  
  184.  
  185.      tonto[9] dap
  186.      DAP> echo=0
  187.      DAP> read /home/miggle/munro/Yorick/include/test1.bas
  188.      DAP> test1(1)
  189.  
  190.     CPU (sec)   SYS (sec)
  191.         8.130        .000
  192.      DAP> test1(10)
  193.  
  194.     CPU (sec)   SYS (sec)
  195.        81.520        .050
  196.      DAP> end
  197.      91.5u 0.6s 3:02 50%
  198.  
  199.  
  200.      tonto[10] yorick
  201.       Yorick ready.  For help type 'help'
  202.      > #include "/home/miggle/munro/Yorick/include/test1.i"
  203.      > test1
  204.             Timing Category     CPU sec  System sec    Wall sec
  205.               Time per pass       0.400       0.000       0.410
  206.              Total time       0.400       0.000       0.410
  207.       -----Total Elapsed Times-----       0.720       0.070      14.810
  208.      > test1,20
  209.             Timing Category     CPU sec  System sec    Wall sec
  210.               Time per pass       0.412       0.000       0.599
  211.              Total time       8.240       0.000      11.980
  212.       -----Total Elapsed Times-----       8.960       0.070      32.350
  213. #endif
  214.  
  215. /* ------------------- Ablation front tracker ---------------------- */
  216.  
  217. /* The ablation front location (rorab) is defined here by the condition that
  218.    the matter temperature reaches a particular fraction of the current drive
  219.    temperature.  Five different hypotheses about this fraction (.4 thru .8)
  220.    are calculated simultaneously.  These can be compared with the ablation
  221.    front location (rorhx) defined by the location of peak phot (X-ray
  222.    absorption), or defined by the steepest gradient (rorsx), or to the
  223.    location of maximum density (rordx).  Other useful parameters are also
  224.    recorded for later perusal.  */
  225. func accum(irec)
  226. {
  227.   /* Input arrays are external -- they would ordinarily be read from
  228.      successive time history dumps of a hydro program.  */
  229.   /* time-independent input quantities */
  230.   extern kmax, lmax;   /* number of points for mesh arrays */
  231.   extern fraction;     /* list of fractions for front definition */
  232.   extern ror;          /* column density coordinate */
  233.  
  234.   /* Since the direction of the ablation is known, the various searches
  235.      begin at the location of the fron on the previous call to accum;
  236.      these are stored in the following two external variables, which
  237.      must be initialized to 37 before the loop calling accum begins.  */
  238.   extern prevScale, prevAbl;
  239.  
  240.   /* time-dependent input quantities */
  241.   extern time;         /* current problem time */
  242.   extern zt;           /* position (point centered mesh array) */
  243.   extern v;            /* velocity (point centered mesh array) */
  244.   extern s;            /* radiation temperature (zone centered mesh array) */
  245.   extern e;            /* matter temperature (zone centered mesh array) */
  246.   extern d;            /* density (zone centered mesh array) */
  247.   extern p;            /* pressure (zone centered mesh array) */
  248.   extern h;            /* heating (zone centered mesh array) */
  249.  
  250.   /* output quantities */
  251.   extern times;        /* list of problem times */
  252.   extern rorab;        /* nfraction-by-ntimes ablation depths (rho*r) */
  253.   extern rordx;        /* ntimes rho*r at peak density */
  254.   extern rorsx;        /* ntimes rho*r at steepest density gradient */
  255.   extern rorhx;        /* ntimes rho*r at peak phot (heating rate) */
  256.   extern dmax;         /* ntimes peak densities */
  257.   extern zs;           /* nfraction-by-ntimes ablation front positions */
  258.   extern vel;          /* ntimes unablated slab velocities */
  259.   extern acc;          /* ntimes unablated slab accelerations */
  260.   extern scl;          /* ntimes density scale lengths at front */
  261.  
  262.   mid = (kmax+1)/2;    /* slab is ablating in l-direction -- calculate
  263.               front for k-zone in middle */
  264.   nfr= numberof(fraction);
  265.   rorzc= ror(zcen);
  266.  
  267.   /* save time */
  268.   times(irec) = time;
  269.  
  270.   /* get choices for ablation temperatures */
  271.   abltemp= fraction*s(mid,lmax);
  272.   epc= e(mid,2:lmax-1);
  273.  
  274.   i= prevAbl-1;  /* epc indices 1 less than e indices */
  275.   for (j=nfr ; j>0 ; j--) {
  276.     ablt= abltemp(j);
  277.     while (epc(i)>ablt) i--;  /* find first zone with lower temperature */
  278.     if (j==nfr) prevAbl= i+1;
  279.     rorab(j,irec)= rorzc(i) +
  280.       (ablt-epc(i))*(rorzc(i+1)-rorzc(i))/(epc(i+1)-epc(i));
  281.   }
  282.  
  283.   /* get peak density and its location */
  284.   dmax(irec)= d(mid, max);
  285.   rordx(irec)= peak(rorzc, d(mid,2:));
  286.  
  287.   /* get location of peak phot heating */
  288.   rorhx(irec) = peak(rorzc(1:35), h(mid,2:36))
  289.  
  290.   /* record position, velocity, and acceleration of zone at peak density */
  291.   i= d(mid, mxx);
  292.   zs(irec)= 0.5*(zt(mid,i)+zt(mid,i-1));
  293.   vel(irec)= 0.5*(v(mid,i)+v(mid,i-1));
  294.   acc(irec)= (0.5*(p(mid,i+1)+p(mid,i-1))-p(mid,i))/(ror(i)-ror(i-1));
  295.  
  296.   /* get density scale height, location of steepest gradient */
  297.   i= min(36, prevScale+4);
  298.   rscl= deriv(d(mid,2:i), rorzc(1:i-1));
  299.   scl(irec)= 1.0/(max(rscl)+1.e-35);
  300.   rorsx(irec) = peak(ror(2:i-1), rscl);
  301.   prevScale= rscl(mxx)+1;
  302.  
  303.   return;
  304. }
  305.  
  306. /* Mock up realistic ablation problem.  */
  307.  
  308. ntimes= 84;
  309. kmax= 11;
  310. lmax= 38;  /* this cannot be changed */
  311.  
  312. time= 1.5;
  313.  
  314. fraction= [.4,.5,.6,.7,.8];  /* MUST be in increasing order */
  315.  
  316. /* the following values represent a very rough fit to an "actual"
  317.    ablating slab problem (that is, the overall shapes of the curves
  318.    are realistic, although the numerical values are not consistent) */
  319. x= span(1,38,38);
  320. y= 4.e-4*atan(x/8) + 1.7e-3*atan((x-16)/3);
  321. y= (2./pi)*(y - y(1));    /* rho*r coordinate */
  322. ror= y(37)-y;             /* ror measured from right (ablating) side */
  323.  
  324. /* v and p were originally fit on a uniform rho*r grid */
  325. yeven= span(min(y), max(y), 38);
  326. tmp= array(0.0, 38);
  327. tmp(1:35)= -0.2 + 0.4*(x(1:35)-x(1))/(x(35)-x(1));
  328. tmp(35:38)= 0.2 + 0.3*(x(35:38)-x(35))/(x(38)-x(35));
  329. v= interp(tmp, yeven, y)(-:1:kmax,);
  330.  
  331. /* original fits to density, pressure, temperature, and heating were for
  332.    point centered data -- the strict uncp algorithm will not work here,
  333.    however, since the input was not actually formed as pairwise averages
  334.    of another array... */
  335. func uncp38(xpc)
  336. {
  337.   x= array(0.0, kmax, 38);
  338.   /* here is "true" uncp algorithm:
  339.      x(2,2)= xpc(1);
  340.      for (i=3 ; i<=38 ; i++) x(2,i)= 2*xpc(i-1)-x(2,i-1);
  341.      x(3:,)= x(2,-,);
  342.    */
  343.   x(2:,2:)= xpc(-,zcen);
  344.   return x;
  345. }
  346.  
  347. dpc= exp(-40./(x+1)^2 - 50./(x-41)^2 + 1.0-0.22*x);
  348. d= uncp38(dpc);
  349.  
  350. tmp= -0.01/(x+7)^2 - 0.0001/(x-41)^2;
  351. tmp= tmp-tmp(1) + (tmp(1)-tmp(38))*(x-x(1))/(x(38)-x(1));
  352. tmp(2:10)= tmp(2) + (tmp(10)-tmp(2))*(x(2:10)-x(2))/(x(10)-x(2));
  353. ppc= interp(tmp, yeven, y);
  354. p= uncp38(ppc);
  355.  
  356. tmp= ppc/dpc;
  357. tmp(19:38)= tmp(19)*(1.0 + 0.1*(x(19:38)-x(19))/(x(38)-x(19)));
  358. tmp(1:5)= max(tmp)/25. + (tmp(5)-max(tmp)/25.)*(x(1:5)-x(1))/(x(5)-x(1));
  359. e= s= uncp38(tmp);
  360.  
  361. tmp(1:9)= 3.6*(x(1:9)-x(1))/(x(9)-x(1));
  362. tmp(9:20)= 3.6 - 2.6*(x(9:20)-x(9))/(x(20)-x(9));
  363. tmp(21:33)= tmp(20);
  364. tmp(33:38)= 1.0 - (x(33:38)-x(33))/(x(38)-x(33));
  365. h= uncp38(tmp);
  366.  
  367. /* declare arrays which will hold results of the calculation */
  368. times= array(0.0, ntimes);
  369. rorab= array(0.0, numberof(fraction), ntimes);
  370. rordx= rorhx= rorsx= acc= vel= scl= zs= dmax= 0.0*times;
  371.  
  372. /* work backwards from d and rhor*r to get zt coordinate */
  373. zt= array(0.0, kmax, 38);
  374. zt(,2:)= (y(dif)/d(2,2:))(-,psum);
  375.  
  376. /* compute parabolic maximum of a curve
  377.      peak = value of xin where yin is maximum   */
  378. func peak(xin,yin)
  379. {
  380.   imax= yin(mxx);
  381.   if (imax==1 || imax==numberof(yin)) return double(imax);
  382.   ymax= yin(imax);
  383.   dy1= ymax-yin(imax-1);
  384.   dy2= yin(imax+1)-ymax;
  385.   xmax= xin(imax);
  386.   dx1= xmax-xin(imax-1);
  387.   dx2= xin(imax+1)-xmax;
  388.   return xmax - (dy2*dx1*dx1+dy1*dx2*dx2)/(dy2*dx1-dy1*dx2);
  389. }
  390.  
  391. /*   deriv(y,x) = derivative of y wrt x
  392.   example:
  393.     acc = deriv( deriv(position,time) , time)
  394.        puts acc equal to the second derivative of position wrt time  */
  395. func deriv(y,x)
  396. {
  397.   return y(dif)/x(dif);
  398. }
  399.